rm(list=ls())

library(tidyverse)
library(marginaleffects)
library(plm)
library(nnet)
library(gridExtra)
library(modelsummary)
library(lme4)


# function to extract legend from plot 
get_only_legend <- function(plot) { 
  plot_table <- ggplot_gtable(ggplot_build(plot)) 
  legend_plot <- which(sapply(plot_table$grobs, function(x) x$name) == "guide-box") 
  legend <- plot_table$grobs[[legend_plot]] 
  return(legend) 
} 

try(setwd("___your directory here___") , silent=TRUE)


pres <- read.csv("data/pres.csv")

pres <- pres |>
  mutate(plotprov=fct_relevel(plotprov,c("Aceh","N. Sumatra","W. Sumatra",
                                         "Riau","Riau Islands","Bangka-Belitung",
                                         "Jambi","S. Sumatra","Bengkulu",
                                         "Lampung","Banten","DKI Jakarta",
                                         "West Java","Central Java","Yogyakarta",
                                         "East Java","Bali","W. Nusa Tenggara",
                                         "E. Nusa Tenggara","W. Kalimantan",
                                         "S. Kalimantan","C. Kalimantan",
                                         "E. Kalimantan","N. Kalimantan",
                                         "N. Sulawesi","Gorontalo","C. Sulawesi",
                                         "W. Sulawesi","S. Sulawesi","S.E. Sulawesi",
                                         "Maluku","N. Maluku","S.W. Papua",
                                         "W. Papua","Papua","C./S./H. Papua")))

  
p <- ggplot(na.omit(pres), aes(x=pg_share, y=ps_share, size=votes, shape=p_winner)) + 
  #geom_smooth(method=lm, se=FALSE, col="gray") + 
  geom_point(alpha=.6) + 
  xlim(0, 1) + ylim(0, 1) + ylab("Prabowo-Sandiaga 2019 Vote Share") +
  xlab("Prabowo-Gibran 2024 Vote Share") +
  labs(size="Total Votes", shape="Prabowo Won\nDistrict (2024)") + 
  geom_hline(yintercept=.5, linetype="dotted") + geom_vline(xintercept=.5, linetype="dotted") +
  scale_size_continuous(range = c(.5,5)) + 
  scale_shape_manual(values=c(1, 16)) + 
  annotate("text", x=.2, y=1, label= "Supported Prabowo in 2019 but not 2024") +
  annotate("text", x=.8, y=1, label= "Supported Prabowo in 2019 and 2024") +
  annotate("text", x=.15, y=0, label= "Opposed Prabowo in 2019\nand lower support in 2024") +
  annotate("text", x=.85, y=0, label= "Opposed Prabowo in 2019\nbut supported Prabowo in 2024") +
  theme_classic() 
p
ggsave(filename = "output/figure1.png", plot = p)
ggsave(filename = "output/figure1.pdf", plot = p)

summary(lm(pg_share ~ ps_share, data=pres))
pres$difference <- pres$pg_share - pres$ps_share

mod1 <- lm(difference ~ muslim + javanese + urban + devindex + ps_share + ethfrac , data=pres)
mod2 <- lm(difference ~ muslim + javanese + urban + devindex + ps_share + first_share , data=pres)
mod3 <- lm(difference ~ muslim + javanese + urban + devindex + ps_share + ethfrac + relfrac , data=pres)

modelsummary(list(mod1,mod2,mod3), stars=TRUE, output="output/appendix_table_A1.docx",
             coef_map=c("muslim","javanese","urban","ethfrac","first_share","relfrac","devindex","ps_share"))

mod5 <- lmer(difference ~ muslim + javanese + urban + devindex + ethfrac + ps_share + 
               (muslim | plotprov), data=pres)
avg_comparisons(mod5, variables = "muslim")

cmp = comparisons(mod5, comparison="dydx",
                  variables = c("muslim"), 
                  newdata = datagrid(plotprov = unique))
p <- ggplot(cmp,
       aes(x = plotprov, y = estimate,
           ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0, linetype = 3) +
  geom_pointrange() + 
  labs(x = "", y = "Effect of Muslim Population Share") +
  theme(axis.text.x = element_text(angle = 90, vjust = .5)) + theme_classic() + 
  theme(axis.text.x = element_text(angle = 90))
ggsave(filename = "output/figureA1.png", plot = p)
ggsave(filename = "output/figureA1.pdf", plot = p)

modelsummary(mod5, stars=TRUE, output="output/appendix_table_A2.docx")



p <- ggplot(na.omit(pres), 
            aes(x=pg_share, y=ps_share, color=muslim, size=votes, shape=p_winner)) +
  geom_point(alpha=.6) + 
  xlim(0, 1) + ylim(0, 1) + ylab("Prabowo-Sandiaga 2019 Vote Share") +
  xlab("Prabowo-Gibran 2024 Vote Share") + 
  labs(size="Total Votes", shape="Prabowo Won\nDistrict (2024)") + labs(color="Proportion Muslim") +
  geom_hline(yintercept=.5, linetype="dotted") + geom_vline(xintercept=.5, linetype="dotted") +
  scale_size_continuous(range = c(.5,5)) + 
  scale_shape_manual(values=c(1, 16)) + 
  theme_classic() + facet_wrap(~plotprov, ncol=6) 
p
ggsave(filename = "output/figure3.png", plot = p)
ggsave(filename = "output/figure3.pdf", plot = p)



islamplot <- pres |>
  group_by(plotprov) |>
  mutate(islamprov = mean(muslim)) %>% 
  ungroup() |>
  filter(islamprov>.9)

p <- ggplot(islamplot, 
            aes(x=ga_share, y=ps_share, color=javanese, size=votes, shape=p_winner)) + 
  geom_point(alpha=.6) + 
  xlim(0, .5) + ylim(0, 1) + ylab("Prabowo-Sandiaga 2019 Vote Share") +
  xlab("Ganjar-Mahfud 2024 Vote Share") + 
  labs(size="Total Votes", shape="Prabowo Won\nDistrict (2024)") + labs(color="Proportion Javanese") +
  geom_hline(yintercept=.5, linetype="dotted") + geom_vline(xintercept=.5, linetype="dotted") +
  scale_size_continuous(range = c(.5,5)) + 
  scale_shape_manual(values=c(1, 16)) + 
  theme_classic() + facet_wrap(~plotprov, ncol=4) 
p
ggsave(filename = "output/figure4.png", plot = p)
ggsave(filename = "output/figure4.pdf", plot = p)



pres <- pres |>
  mutate(third = pmin(Ganjar,Anies,Prabowo),
         second = if_else(Anies > Ganjar & Ganjar > Prabowo, Ganjar,
                                ifelse(Anies > Prabowo & Prabowo > Ganjar, Prabowo,
                                       ifelse(Ganjar > Prabowo & Prabowo > Anies,Prabowo,
                                              ifelse(Ganjar > Anies & Anies > Prabowo,Anies,
                                                     ifelse(Prabowo > Anies & Anies > Ganjar,Anies,
                                                            ifelse(Prabowo > Ganjar & Ganjar > Anies,Ganjar,NA)))))),
         loser = case_when(
           Anies < Ganjar & Anies < Prabowo ~ "Anies",
           Prabowo < Ganjar & Prabowo < Anies ~ "Prabowo",
           Ganjar < Prabowo & Ganjar < Anies ~ "Ganjar"),
         lasttwo = third/second)


p <- ggplot(pres[!is.na(pres$loser),], aes(lasttwo)) + 
  geom_histogram(binwidth = .1) + 
  ylab("Number of Districts") + xlab("Ratio of Third- to Second-Place Votes") +
  theme_classic() + facet_wrap(~loser)
ggsave(filename = "output/figure5.png", plot = p, width=6.5, height=3)
ggsave(filename = "output/figure5.pdf", plot = p, width=6.5, height=3)
